home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / pl3d.i < prev    next >
Text File  |  1996-02-13  |  28KB  |  906 lines
  1. /*
  2.    PL3D.I
  3.    Viewing transforms and other aids for 3D plotting.
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1996.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. /* General overview:
  11.  
  12.    (1) Viewing transform machinery.  Arguably the simplest model
  13.        is the CAD/CAM notion that the object you see is oriented
  14.        as you see it in the current picture.  You can then move
  15.        it left, right, up, down, or toward or away from you,
  16.        or you can rotate it about any of the three axes (horizontal,
  17.        vertical, or out of the screen).  The xyz coordinates of the
  18.        object remains unchanged throughout all of this, but this
  19.        object coordinate system changes relative to the fixed
  20.        xyz of the viewer, in which x is always to the right, y is
  21.        up, and z is directed out of the screen.  Initially, the
  22.        two coordinate systems coincide.
  23.  
  24.        rot3, xangle,yangle,zangle
  25.          Rotate the object about viewer's x-axis by xangle, then
  26.      about viewer's y-axis by yangle, then about viewer's
  27.      z-axis by zangle
  28.        mov3, xchange,ychange,zchange
  29.          Move the object by the specified amounts.
  30.  
  31.        setz3, zcamera
  32.          The "camera" is located at (0,0,zcamera) in the viewer's
  33.      coordinate system, looking in the minus-z direction.
  34.      Initially, zcamera is very large, and the magnification
  35.      factor is correspondingly large, giving an isometric view.
  36.      Decreasing zcamera makes the perspective more extreme.
  37.      If parts of the object are behind the camera, strange things
  38.      may happen.
  39.  
  40.        undo3
  41.        undo3, n
  42.          Undo the last N (default 1) viewpoint commands (rot3, mov3,
  43.      or setz3).  Up to 100 viewpoint changes are remembered.
  44.        viewpoint= save3()
  45.        ...
  46.        restore3, viewpoint
  47.          The current viewpoint transformation can be saved and later
  48.      restored.
  49.  
  50.        gnomon, on_off
  51.          Toggle the gnomon (a simple display showing the orientation
  52.      of the xyz axes of the object).
  53.  */
  54.  
  55. /* ------------------------------------------------------------------------ */
  56.  
  57. func rot3(xa,ya,za)
  58. /* DOCUMENT rot3, xa,ya,za
  59.      rotate the current 3D plot by XA about viewer's x-axis,
  60.      YA about viewer's y-axis, and ZA about viewer's z-axis.
  61.    SEE ALSO: orient3, mov3, aim3, setz3, undo3, save3, restore3, light3
  62.  */
  63. {
  64.   if (is_void(xa)) xa= 0.;
  65.   if (is_void(ya)) ya= 0.;
  66.   if (is_void(za)) za= 0.;
  67.   x= [1.,0.,0.];
  68.   y= [0.,1.,0.];
  69.   z= [0.,0.,1.];
  70.   _rot3, za, x, y;
  71.   _rot3, ya, z, x;
  72.   _rot3, xa, y, z;
  73.   _setrot3, [x,y,z](,+)*_getrot3()(+,);
  74. }
  75.  
  76. func _rot3(a, &x, &y)
  77. {
  78.   ca= cos(a);
  79.   sa= sin(a);
  80.   a= x;
  81.   x=  ca*a + sa*y;
  82.   y= -sa*a + ca*y;
  83. }
  84.  
  85. func mov3(xa,ya,za)
  86. /* DOCUMENT mov3, xa,ya,za
  87.      move the current 3D plot by XA along viewer's x-axis,
  88.      YA along viewer's y-axis, and ZA along viewer's z-axis.
  89.    SEE ALSO: rot3, orient3, setz3, undo3, save3, restore3, light3
  90.  */
  91. {
  92.   if (is_void(xa)) xa= 0.;
  93.   if (is_void(ya)) ya= 0.;
  94.   if (is_void(za)) za= 0.;
  95.   _setorg3, _getorg3() - _getrot3()(+,)*[xa,ya,za](+);
  96. }
  97.  
  98. func aim3(xa,ya,za)
  99. /* DOCUMENT aim3, xa,ya,za
  100.      move the current 3D plot to put the point (XA,YA,ZA) in object
  101.      coordinates at the point (0,0,0) -- the aim point -- in the
  102.      viewer's coordinates.  If any of XA, YA, or ZA is nil, it defaults
  103.      to zero.
  104.    SEE ALSO: mov3, rot3, orient3, setz3, undo3, save3, restore3, light3
  105.  */
  106. {
  107.   if (is_void(xa)) xa= 0.;
  108.   if (is_void(ya)) ya= 0.;
  109.   if (is_void(za)) za= 0.;
  110.   _setorg3, [xa,ya,za];
  111. }
  112.  
  113. func setz3(zc)
  114. /* DOCUMENT setz3, zc
  115.      Set the camera position to z=ZC (x=y=0) in the viewer's coordinate
  116.      system.  If ZC is nil, set the camera to infinity (default).
  117.    SEE ALSO: rot3, orient3, undo3, save3, restore3, light3
  118.  */
  119. {
  120.   if (!is_void(zc)) {
  121.     if (dimsof(zc)(1)) error, "camera position must be scalar";
  122.     zc= double(zc);
  123.   }
  124.   _setzc3, zc;
  125. }
  126.  
  127. func orient3(phi, theta)
  128. /* DOCUMENT orient3, phi, theta
  129.          or orient3, phi
  130.          or orient3, , theta
  131.          or orient3
  132.      Set the "orientation" of the object to (PHI,THETA).  "Orientations"
  133.      are a subset of the possible rotation matrices in which the z-axis
  134.      of the object appears vertical on the screen (that is, the object
  135.      z-axis projects onto the viewer y-axis).  The THETA angle is the
  136.      angle from the viewer y-axis to the object z-axis, positive if
  137.      the object z-axis is tilted toward you (toward viewer +z).  PHI is
  138.      zero when the object x-axis coincides with the viewer x-axis.  If
  139.      neither PHI nor THETA is specified, PHI defaults to -pi/4 and
  140.      THETA defaults to pi/6.  If only one of PHI or THETA is specified,
  141.      the other remains unchanged, unless the current THETA is near pi/2,
  142.      in which case THETA returns to pi/6, or unless the current
  143.      orientation does not have a vertical z-axis, in which case the
  144.      unspecified value returns to its default.
  145.  
  146.      Unlike rot3, orient3 is not a cumulative operation.
  147.  
  148.    SEE ALSO: rot3, mov3, aim3, save3, restore3, light3
  149.  */
  150. {
  151.   if (is_void(theta) && is_void(phi)) {
  152.     theta= _orient3_theta;
  153.     phi= _orient3_phi;
  154.   } else if (is_void(theta) || is_void(phi)) {
  155.     z= _getrot3()(,+)*[0.,0.,1.](+);
  156.     if (abs(z(1))>1.e-6) {
  157.       /* object z-axis not aligned with viewer y-axis */
  158.       if (is_void(theta)) theta= _orient3_theta;
  159.       else phi= _orient3_phi;
  160.     } else if (is_void(theta)) {
  161.       if (abs(z(2))<1.e-6) theta= _orient3_theta;
  162.       else theta= atan(z(3),z(2));
  163.     } else /*if (is_void(phi))*/ {
  164.       y= [0.,z(3),-z(2)];  /* in object xy-plane */
  165.       x= _getrot3()(,+)*[1.,0.,0.](+);
  166.       phi= atan(y(+)*x(+), x(1));
  167.     }
  168.   }
  169.   x= [1.,0.,0.];
  170.   y= [0.,1.,0.];
  171.   z= [0.,0.,1.];
  172.   _rot3, theta, y, z;
  173.   _rot3, phi, z, x;
  174.   _setrot3, [x,-z,y];
  175. }
  176.  
  177. /* unless user has supplied alternative defaults, set orient3 defaults */
  178. if (is_void(_orient3_theta)) _orient3_theta= pi/6;
  179. if (is_void(_orient3_phi)) _orient3_phi= -pi/4;
  180.  
  181. func save3(void)
  182. /* DOCUMENT view= save3()
  183.      Save the current 3D viewing transformation and lighting.
  184.    SEE ALSO: restore3, rot3, mov3, aim3, light3
  185.  */
  186. {
  187.   return _cpy(_draw3_list, _draw3_n);
  188. }
  189.  
  190. func restore3(view)
  191. /* DOCUMENT restore3, view
  192.      Restore a previously saved 3D viewing transformation and lighting.
  193.      If VIEW is nil, rotate object to viewer's coordinate system.
  194.    SEE ALSO: restore3, rot3, mov3, aim3, light3
  195.  */
  196. {
  197.   if (!is_void(view)) view= _cpy(view);
  198.   else view= _cat(_cpy(_draw3_view), _cpy(_light3_list));
  199.   _cdr, view, _draw3_n, _cdr(_draw3_list, _draw3_n, []);
  200.   _draw3_list= view;
  201.   _undo3_set, restore3, old;
  202. }
  203.  
  204. /* set default viewing direction if user hasn't already done so */
  205. if (is_void(_draw3_view)) _draw3_view= _lst(unit(3), [0.,0.,0.], []);
  206. _draw3_nv= _len(_draw3_view);
  207.  
  208. /* ------------------------------------------------------------------------ */
  209.  
  210. func light3(ambient=,diffuse=,specular=,spower=,sdir=)
  211. /* DOCUMENT light3, ambient=a_level,
  212.                     diffuse=d_level,
  213.             specular=s_level,
  214.             spower=n,
  215.             sdir=xyz
  216.      Sets lighting properties for 3D shading effects.
  217.      A surface will be shaded according to its to its orientation
  218.      relative to the viewing direction.
  219.  
  220.      The ambient level A_LEVEL is a light level (arbitrary units)
  221.      that is added to every surface independent of its orientation.
  222.  
  223.      The diffuse level D_LEVEL is a light level which is proportional
  224.      to cos(theta), where theta is the angle between the surface
  225.      normal and the viewing direction, so that surfaces directly
  226.      facing the viewer are bright, while surfaces viewed edge on are
  227.      unlit (and surfaces facing away, if drawn, are shaded as if they
  228.      faced the viewer).
  229.  
  230.      The specular level S_LEVEL is a light level proportional to a high
  231.      power spower=N of 1+cos(alpha), where alpha is the angle between
  232.      the specular reflection angle and the viewing direction.  The light
  233.      source for the calculation of alpha lies in the direction XYZ (a
  234.      3 element vector) in the viewer's coordinate system at infinite
  235.      distance.  You can have ns light sources by making S_LEVEL, N, and
  236.      XYZ (or any combination) be vectors of length ns (3-by-ns in the
  237.      case of XYZ).  (See source code for specular_hook function
  238.      definition if powers of 1+cos(alpha) aren't good enough for you.)
  239.  
  240.      With no arguments, return to the default lighting.
  241.  
  242.    EXAMPLES:
  243.      light3, diffuse=.1, specular=1., sdir=[0,0,-1]
  244.        (dramatic "tail lighting" effect)
  245.      light3, diffuse=.5, specular=1., sdir=[1,.5,1]
  246.        (classic "over your right shoulder" lighting)
  247.      light3, ambient=.1,diffuse=.1,specular=1.,
  248.              sdir=[[0,0,-1],[1,.5,1]],spower=[4,2]
  249.        (two light sources combining previous effects)
  250.  
  251.    SEE ALSO: rot3, save3, restore3
  252.  */
  253. {
  254.   old= _cpy(_cdr(_draw3_list,_draw3_nv),5);
  255.  
  256.   flags= 0;
  257.   if (!is_void(ambient)) {
  258.     if (dimsof(ambient)(1)) error, "ambient light level must be scalar";
  259.     flags|= 1;
  260.     _car, _draw3_list, _draw3_nv+1, double(ambient);
  261.   }
  262.   if (!is_void(diffuse)) {
  263.     if (dimsof(diffuse)(1)) error, "diffuse light level must be scalar";
  264.     flags|= 2;
  265.     _car, _draw3_list, _draw3_nv+2, double(diffuse);
  266.   }
  267.  
  268.   if (!is_void(specular)) flags|= 4;
  269.   else specular= _car(_draw3_list, _draw3_nv+3);
  270.   if (!is_void(spower)) flags|= 8;
  271.   else spower= _car(_draw3_list, _draw3_nv+4);
  272.   if (!is_void(sdir)) {
  273.     dims= dimsof(sdir);
  274.     if (dims(1)<1 || dims(2)!=3)
  275.       error, "lighting direction must be 3 vector or 3-by-ns array"
  276.     flags|= 16;
  277.   } else {
  278.     sdir= _car(_draw3_list, _draw3_nv+5);
  279.   }
  280.   if (flags&28) {
  281.     if (is_void(dimsof(specular,spower,sdir(1,..))))
  282.       error, "specular, spower, and sdir not conformable";
  283.     if (flags&4) _car, _draw3_list, _draw3_nv+3, double(specular);
  284.     if (flags&8) _car, _draw3_list, _draw3_nv+4, spower;
  285.     if (flags&16) _car, _draw3_list, _draw3_nv+5, double(sdir);
  286.   }
  287.  
  288.   if (!flags) {
  289.     for (i=1 ; i<=5 ; ++i)
  290.       _car, _draw3_list, _draw3_nv+i, _car(_light3_list,i);
  291.   }
  292.  
  293.   _undo3_set, _light3, old;
  294. }
  295.  
  296. func _light3(arg)
  297. {
  298.   for (i=1 ; i<=5 ; ++i)
  299.     _car, _draw3_list, _draw3_nv+i, _car(arg,i);
  300. }
  301.  
  302. /* set default values if user hasn't already done so */
  303. if (is_void(_light3_ambient)) _light3_ambient= 0.2;
  304. if (is_void(_light3_diffuse)) _light3_diffuse= 1.0;
  305. if (is_void(_light3_specular)) _light3_specular= 0.0;
  306. if (is_void(_light3_spower)) _light3_spower= 2;
  307. if (is_void(_light3_sdir)) _light3_sdir= [1.0, 0.5, 1.0]/sqrt(2.25);
  308. _light3_list= _lst(_light3_ambient,_light3_diffuse,_light3_specular,
  309.            _light3_spower,_light3_sdir);
  310.  
  311. func get3_light(xyz, nxyz)
  312. /* DOCUMENT get3_light(xyz, nxyz)
  313.          or get3_light(xyz)
  314.  
  315.      return 3D lighting for polygons with vertices XYZ.  If NXYZ is
  316.      specified, XYZ should be 3-by-sum(nxyz), with NXYZ being the
  317.      list of numbers of vertices for each polygon (as for the plfp
  318.      function).  If NXYZ is not specified, XYZ should be a quadrilateral
  319.      mesh, 3-by-ni-by-nj (as for the plf function).  In the first case,
  320.      the return value is numberof(NXYZ); in the second case, the
  321.      return value is (ni-1)-by-(nj-1).
  322.  
  323.      The parameters of the lighting calculation are set by the
  324.      light3 function.
  325.  
  326.    SEE ALSO: light3, set3_object, get3_normal, get3_centroid
  327.  */
  328. {
  329.   list= _cdr(_draw3_list, _draw3_nv);
  330.   ambient= _nxt(list);
  331.   diffuse= _nxt(list);
  332.   specular= _nxt(list);
  333.   spower= _nxt(list);
  334.   sdir= _nxt(list);
  335.  
  336.   /* get normal */
  337.   normal= get3_normal(xyz, nxyz);
  338.  
  339.   /* get direction to viewer's eye (camera) */
  340.   zc= _getzc3();
  341.   if (is_void(zc)) {
  342.     view= [0.,0.,1.];
  343.   } else {
  344.     view= [0.,0.,zc]-get3_centroid(xyz, nxyz);
  345.     m1= abs(view(1,..),view(2,..),view(3,..))(-,..);
  346.     m1= m1 + (m1==0.0);
  347.     view/= m1;
  348.   }
  349.  
  350.   /* do lighting calculation */
  351.   nv= (normal*view)(sum,..);
  352.   light= ambient + diffuse*abs(nv);
  353.   if (anyof(specular)) {
  354.     sdir= sdir(,*);
  355.     sdir/= abs(sdir(1,),sdir(2,),sdir(3,))(-,);
  356.     sv= sdir(+,*)*view(+,..);
  357.     sn= sdir(+,*)*normal(+,..);
  358.     m1= max(sn*nv(-,) - 0.5*sv + 0.5, 1.e-30);  /* max(1+cos(alpha),0) */
  359.     if (is_func(specular_hook))
  360.       m1= specular_hook(m1, abs(nv)(-,..), spower);
  361.     else
  362.       m1= m1^spower;
  363.     light+= (specular(*)*m1)(sum,..);
  364.   }
  365.  
  366.   return light;
  367. }
  368.  
  369. func get3_normal(xyz, nxyz)
  370. /* DOCUMENT get3_normal(xyz, nxyz)
  371.          or get3_normal(xyz)
  372.  
  373.      return 3D normals for polygons with vertices XYZ.  If NXYZ is
  374.      specified, XYZ should be 3-by-sum(nxyz), with NXYZ being the
  375.      list of numbers of vertices for each polygon (as for the plfp
  376.      function).  If NXYZ is not specified, XYZ should be a quadrilateral
  377.      mesh, 3-by-ni-by-nj (as for the plf function).  In the first case,
  378.      the return value is 3-by-numberof(NXYZ); in the second case, the
  379.      return value is 3-by-(ni-1)-by-(nj-1).
  380.  
  381.      The normals are constructed from the cross product of the lines
  382.      joining the midpoints of two edges which as nearly quarter the
  383.      polygon as possible (the medians for a quadrilateral).  No check
  384.      is made that these not be parallel; the returned "normal" is
  385.      [0,0,0] in that case.  Also, if the polygon vertices are not
  386.      coplanar, the "normal" has no precisely definable meaning.
  387.  
  388.    SEE ALSO: get3_centroid, get3_light
  389.  */
  390. {
  391.   if (is_void(nxyz)) {
  392.     /* if no polygon list is given, assume xyz is 2D mesh */
  393.     /* form normal as cross product of medians */
  394.     m1= xyz(,zcen,dif);
  395.     m2= xyz(,dif,zcen);
  396.  
  397.   } else {
  398.     /* with polygon list, more elaborate calculation required */
  399.     frst= nxyz(psum)-nxyz+1;
  400.  
  401.     /* form normal by getting two approximate diameters
  402.      * (reduces to above medians for quads) */
  403.     n2= (nxyz+1)/2;
  404.     zero= frst-1;
  405.     c0= 0.5*(xyz(,zero+1)+xyz(,zero+2));
  406.     i= zero+n2;  /* n2>=2, nxyz>=3 */
  407.     c1= 0.5*(xyz(,i)+xyz(,i+1));
  408.     i= 1+n2/2;
  409.     c2= 0.5*(xyz(,zero+i)+xyz(,zero+i+1));
  410.     i= min(i+n2, nxyz);
  411.     c3= 0.5*(xyz(,zero+i)+xyz(,zero+i%nxyz+1));
  412.     m1= c1 - c0;
  413.     m2= c3 - c2;
  414.   }
  415.  
  416.   /* poly normal is cross product of two medians (or diameters) */
  417.   normal= m1;
  418.   normal(1,..)= n1= m1(2,..)*m2(3,..) - m1(3,..)*m2(2,..);
  419.   normal(2,..)= n2= m1(3,..)*m2(1,..) - m1(1,..)*m2(3,..);
  420.   normal(3,..)= n3= m1(1,..)*m2(2,..) - m1(2,..)*m2(1,..);
  421.   m1= abs(n1,n2,n3)(-,..);
  422.   m1= m1 + (m1==0.0);
  423.   normal/= m1;
  424.  
  425.   return normal;
  426. }
  427.  
  428. func get3_centroid(xyz, nxyz)
  429. /* DOCUMENT get3_centroid(xyz, nxyz)
  430.          or get3_centroid(xyz)
  431.  
  432.      return 3D centroids for polygons with vertices XYZ.  If NXYZ is
  433.      specified, XYZ should be 3-by-sum(nxyz), with NXYZ being the
  434.      list of numbers of vertices for each polygon (as for the plfp
  435.      function).  If NXYZ is not specified, XYZ should be a quadrilateral
  436.      mesh, 3-by-ni-by-nj (as for the plf function).  In the first case,
  437.      the return value is 3-by-numberof(NXYZ); in the second case, the
  438.      return value is 3-by-(ni-1)-by-(nj-1).
  439.  
  440.      The centroids are constructed as the mean value of all vertices
  441.      of each polygon.
  442.  
  443.    SEE ALSO: get3_normal, get3_light
  444.  */
  445. {
  446.   if (is_void(nxyz)) {
  447.     /* if no polygon list is given, assume xyz is 2D mesh */
  448.     centroid= xyz(,zcen,zcen);
  449.  
  450.   } else {
  451.     /* with polygon list, more elaborate calculation required */
  452.     last= nxyz(psum);
  453.     list= histogram(1+last)(1:-1);
  454.     list(1)+= 1;
  455.     list= list(psum);
  456.     centroid= array(0.0, 3, numberof(nxyz));
  457.     centroid(1,)= histogram(list, xyz(1,));
  458.     centroid(2,)= histogram(list, xyz(2,));
  459.     centroid(3,)= histogram(list, xyz(3,));
  460.     centroid/= double(nxyz);
  461.   }
  462.  
  463.   return centroid;
  464. }
  465.  
  466. /* ------------------------------------------------------------------------ */
  467.  
  468. func get3_xy(xyz, &x, &y, &z, getz)
  469. /* DOCUMENT get3_xy, xyz, x, y
  470.          or get3_xy, xyz, x, y, z, 1
  471.  
  472.      Given 3-by-anything coordinates XYZ, return X and Y in viewer's
  473.      coordinate system (set by rot3, mov3, orient3, etc.).  If the
  474.      fifth argument is present and non-zero, also return Z (for use
  475.      in sort3d or get3_light, for example).  If the camera position
  476.      has been set to a finite distance with setz3, the returned
  477.      coordinates will be tangents of angles for a perspective
  478.      drawing (and Z will be scaled by 1/zc).
  479.  
  480.    SEE ALSO: sort3d, get3_light, rot3, setz3, set3_object
  481.  */
  482. {
  483.   /* rotate and translate to viewer's coordinate system */
  484.   xyz= _getrot3()(,+)*(xyz - _getorg3())(+,..);
  485.   x= xyz(1,..);
  486.   y= xyz(2,..);
  487.  
  488.   /* do optional perspective projection */
  489.   zc= _getzc3();
  490.   if (!is_void(zc)) {
  491.     z= xyz(3,..);
  492.     zc= max(zc-z, 0.0);  /* protect behind camera */
  493.     zc+= (zc==0.0)*1.e-35;       /* avoid zero divide */
  494.     x/= zc;
  495.     y/= zc;
  496.     if (getz) z/= zc;
  497.   } else if (getz) {
  498.     z= xyz(3,..);
  499.   }
  500. }
  501.  
  502. func _getrot3(void)
  503. {
  504.   return _car(_draw3_list, 1);
  505. }
  506. func _getorg3(void)
  507. {
  508.   return _car(_draw3_list, 2);
  509. }
  510. func _getzc3(void)
  511. {
  512.   return _car(_draw3_list, 3);
  513. }
  514.  
  515. func _setrot3(x)
  516. {
  517.   _undo3_set, _setrot3, _car(_draw3_list, 1, x);
  518. }
  519. func _setorg3(x)
  520. {
  521.   _undo3_set, _setorg3, _car(_draw3_list, 2, x);
  522. }
  523. func _setzc3(x)
  524. {
  525.   _undo3_set, _setzc3, _car(_draw3_list, 3, x);
  526. }
  527.  
  528. func _undo3_set(fnc, arg)
  529. {
  530.   if (!_in_undo3) {
  531.     if (_len(_undo3_list)>=2*_undo3_limit)
  532.       _cdr, _undo3_list, 2*_undo3_limit-2, [];
  533.     _undo3_list= _cat(_lst(fnc,arg), _undo3_list);
  534.   }
  535.   draw3_trigger;
  536. }
  537.  
  538. _undo3_limit= 100;
  539.  
  540. func undo3(n)
  541. /* DOCUMENT undo3
  542.          or undo3, n
  543.      Undo the effects of the last N (default 1) rot3, orient3, mov3, aim3,
  544.      setz3, or light3 commands.
  545.  */
  546. {
  547.   if (is_void(n)) n= 1;
  548.   n= 2*(n-1);
  549.   if (n<0 || n>_len(_undo3_list))
  550.     error, "not that many items in undo list";
  551.   _in_undo3= 1;  /* flag to skip _undo3_set */
  552.   /* perhaps should save discarded items in a redo list? */
  553.   if (n) _undo3_list= _cdr(_undo3_list, n);
  554.   for (; n>=0 ; n-=2) {
  555.     fnc= _nxt(_undo3_list);
  556.     arg= _nxt(_undo3_list);
  557.     fnc, arg;
  558.   }
  559.   draw3_trigger;
  560. }
  561.  
  562. func set3_object(fnc, arg)
  563. /* DOCUMENT set3_object, drawing_function, _lst(arg1,arg2,...)
  564.  
  565.      set up to trigger a call to draw3, adding a call to the
  566.      3D display list of the form:
  567.  
  568.         DRAWING_FUNCTION, _lst(ARG1, ARG2, ...)
  569.  
  570.      When draw3 calls DRAWING_FUNCTION, the external variable _draw3
  571.      will be non-zero, so DRAWING_FUNCTION can be written like this:
  572.  
  573.      func drawing_function(arg1,arg2,...)
  574.      {
  575.        require, "pl3d.i";
  576.        if (_draw3) {
  577.          list= arg1;
  578.      arg1= _nxt(list);
  579.      arg2= _nxt(list);
  580.      ...
  581.      ...<calls to get3_xy, sort3d, get3_light, etc.>...
  582.      ...<calls to graphics functions plfp, plf, etc.>...
  583.      return;
  584.        }
  585.        ...<verify args>...
  586.        ...<do orientation and lighting independent calcs>...
  587.        set3_object, drawing_function, _lst(arg1,arg2,...);
  588.      }
  589.  
  590.   SEE ALSO: get3_xy, get3_light, sort3d
  591.  */
  592. {
  593.   _draw3_list= _cat(_draw3_list, _lst(fnc,arg));
  594.   draw3_trigger;
  595. }
  596.  
  597. /* ------------------------------------------------------------------------ */
  598.  
  599. func draw3(called_as_idler)
  600. /* DOCUMENT draw3
  601.      Draw the current 3D display list.
  602.      (Ordinarily triggered automatically when the drawing changes.)
  603.  */
  604. {
  605.   if (_draw3_changes) {
  606.     if (called_as_idler) fma;
  607.     /* the first _draw3_n elements of _draw3_list are the viewing
  608.      * transforms, lighting, etc.
  609.      * thereafter, elements are (function,argument-list) pairs
  610.      * the _draw3 flag alerts the functions that these are the draw
  611.      * calls rather than the interactive setup calls */
  612.     limits, square=1;
  613.     _draw3= 1;
  614.     for (list=_cdr(_draw3_list, _draw3_n) ; list ; list=_cdr(list)) {
  615.       fnc= _car(list);
  616.       list= _cdr(list);
  617.       fnc, _car(list);
  618.     }
  619.     if (_gnomon) _gnomon_draw;
  620.     _draw3_changes= [];
  621.   }
  622. }
  623.  
  624. _draw3_list= _cat(_cpy(_draw3_view), _cpy(_light3_list));
  625. _draw3_n= _len(_draw3_list);
  626.  
  627. func draw3_trigger
  628. {
  629.   /* arrange to call draw3 when everything else is finished */
  630.   set_idler, _draw3_idler;
  631.   extern _draw3_changes;
  632.   _draw3_changes= 1;
  633. }
  634.  
  635. func _draw3_idler
  636. {
  637.   draw3, 1;
  638. }
  639.  
  640. func clear3(void)
  641. /* DOCUMENT clear3
  642.      Clear the current 3D display list.
  643.  */
  644. {
  645.   _cdr, _draw3_list, _draw3_n, [];
  646. }
  647.  
  648. func window3(n)
  649. /* DOCUMENT window3
  650.          or window3, n
  651.      initialize style="nobox.gs" window for 3D graphics
  652.  */
  653. {
  654.   window, n, wait=1, style="nobox.gs", legends=0;
  655. }
  656.  
  657. /* ------------------------------------------------------------------------ */
  658.  
  659. func gnomon(on)
  660. /* DOCUMENT gnomon
  661.          or gnomon, onoff
  662.      Toggle the gnomon display.  If ONOFF is non-nil and non-zero,
  663.      turn on the gnomon.  If ONOFF is zero, turn off the gnomon.
  664.  
  665.      The gnomon shows the X, Y, and Z axis directions in the
  666.      object coordinate system.  The directions are labeled.
  667.      The gnomon is always infinitely far behind the object
  668.      (away from the camera).
  669.  
  670.      There is a mirror-through-the-screen-plane ambiguity in the
  671.      display which is resolved in two ways: (1) The (X,Y,Z)
  672.      coordinate system is right-handed, and (2) If the tip of an
  673.      axis projects into the screen, it's label is drawn in opposite
  674.      polarity to the other text on the screen.
  675.  */
  676. {
  677.   old= _gnomon;
  678.   if (is_void(on)) _gnomon~= 1;
  679.   else if (on) _gnomon= 1;
  680.   else _gnomon= 0;
  681.   if (old!=_gnomon) draw3_trigger;
  682. }
  683.  
  684. if (is_void(_gnomon)) _gnomon= 0;
  685.  
  686. func _gnomon_draw(void)
  687. {
  688.   o= [0.,0.,0.];
  689.   x1= [1.,0.,0.];
  690.   y1= [0.,1.,0.];
  691.   z1= [0.,0.,1.];
  692.   xyz= _getrot3()(,+)*[[o,x1],[o,y1],[o,z1]](+,,,);
  693.   xyz*= 0.0013*_gnomon_scale;
  694.   x1= xyz(1,,);
  695.   y1= xyz(2,,);
  696.   z1= xyz(3,2,);
  697.   x0= x1(1,);
  698.   x1= x1(2,);
  699.   y0= y1(1,);
  700.   y1= y1(2,);
  701.   wid= min(_gnomon_scale/18.,6.);
  702.   if (wid<0.5) wid= 0.0;
  703.   plsys, 0;
  704.   pldj, x0+_gnomon_x, y0+_gnomon_y, x1+_gnomon_x, y1+_gnomon_y,
  705.     width=wid, type=1, legend=string(0);
  706.   plsys, 1;
  707.  
  708.   /* compute point size of labels (1/3 of axis length) */
  709.   pts= [8,10,12,14,18,24](digitize(_gnomon_scale/3.0,
  710.                    [9,11,13,16,21]));
  711.   if (_gnomon_scale < 21.) {
  712.     x1*= 21./_gnomon_scale;
  713.     y1*= 21./_gnomon_scale;
  714.   }
  715.   /* label positions: first find shortest axis */
  716.   xy= abs(x1,y1);
  717.   i= xy(mnx);
  718.   jk= [[2,3],[3,1],[1,2]](,i);
  719.   if (xy(i)<1.e-7*xy(sum)) {  /* guarantee not exactly zero */
  720.     x1(i)= -1.e-6*x1(jk)(sum);
  721.     y1(i)= -1.e-6*y1(jk)(sum);
  722.     xy(i)= abs(x1(i),y1(i));
  723.   }
  724.   xyi= xy(i);
  725.   /* next find axis nearest to shortest */
  726.   j= jk(1);
  727.   k= jk(2);
  728.   if (abs(x1(j)*y1(i)-y1(j)*x1(i))*xy(k) >
  729.       abs(x1(k)*y1(i)-y1(k)*x1(i))*xy(j)) {
  730.     jk= j;  j= k;  k= jk;
  731.   }
  732.   /* furthest axis first - move perpendicular to nearer axis */
  733.   xk= -y1(j);
  734.   yk= x1(j);
  735.   xy= abs(xk,yk);
  736.   xk/= xy;
  737.   yk/= xy;
  738.   if (xk*x1(k)+yk*y1(k) < 0.0) { xk= -xk;  yk= -yk; }
  739.   /* nearer axis next - move perpendicular to furthest axis */
  740.   xj= -y1(k);
  741.   yj= x1(k);
  742.   xy= abs(xj,yj);
  743.   xj/= xy;
  744.   yj/= xy;
  745.   if (xj*x1(j)+yj*y1(j) < 0.0) { xj= -xj;  yj= -yj; }
  746.   /* shortest axis last - move perpendicular to nearer */
  747.   xi= -y1(j);
  748.   yi= x1(j);
  749.   xy= abs(xi,yi);
  750.   xi/= xy;
  751.   yi/= xy;
  752.   if (xi*x1(i)+yi*y1(i) < 0.0) { xi= -xi;  yi= -yi; }
  753.  
  754.   /* shortest axis label may need adjustment */
  755.   d= 0.0013*pts;
  756.   if (xyi < d) {
  757.     /* just center it in correct quadrant */
  758.     jk= sign(xi*xj+yi*yj);
  759.     yi= sign(xi*xk+yi*yk);
  760.     xi= jk*xj + yi*xk;
  761.     yi= jk*yj + yi*yk;
  762.     jk= abs(xi, yi);
  763.     xi/= jk;
  764.     yi/= jk;
  765.   }
  766.  
  767.   x= y= [0.,0.,0.];
  768.   x([i,j,k])= [xi,xj,xk];
  769.   y([i,j,k])= [yi,yj,yk];
  770.   x*= d;
  771.   y*= d;
  772.   x+= x1 + _gnomon_x;
  773.   y+= y1 + _gnomon_y;
  774.   chr= ["X","Y","Z"];
  775.   _gnomon_text, chr(i), x(i),y(i), pts, z1(i)<-1.e-6;
  776.   _gnomon_text, chr(j), x(j),y(j), pts, z1(j)<-1.e-6;
  777.   _gnomon_text, chr(k), x(k),y(k), pts, z1(k)<-1.e-6;
  778. }
  779.  
  780. /* recommended _gnomon_scale: 24, 30, 36, 42, 54, or 72 */
  781. if (is_void(_gnomon_scale))
  782.   _gnomon_scale= 30.   /* X,Y,Z axis lengths in points */
  783. if (is_void(_gnomon_x))
  784.   _gnomon_x= 0.18;     /* gnomon origin in system 0 coordinates */
  785. if (is_void(_gnomon_y))
  786.   _gnomon_y= 0.42;
  787.  
  788. func _gnomon_text(chr, x, y, pts, invert)
  789. {
  790.   /* pts= 8, 10, 12, 14, 18, or 24 */
  791.   col= "fg";
  792.   if (invert) {
  793.     plsys, 0;
  794.     plg,[y,y],[x,x],
  795.       type=1,width=2.2*pts,color=col,marks=0,legend=string(0);
  796.     plsys, 1;
  797.     col= "bg";
  798.   }
  799.   plt, chr, x,y,
  800.     justify="CH",color=col,height=pts,font="helvetica",opaque=0;
  801. }
  802.  
  803. /* ------------------------------------------------------------------------ */
  804.  
  805. func sort3d(z, npolys, &list, &vlist)
  806. /* DOCUMENT sort3d(z, npolys, &list, &vlist)
  807.  
  808.      given Z and NPOLYS, with numberof(Z)==sum(npolys), return
  809.      LIST and VLIST such that Z(VLIST) and NPOLYS(LIST) are
  810.      sorted from smallest average Z to largest average Z, where
  811.      the averages are taken over the clusters of length NPOLYS.
  812.      Within each cluster (polygon), the cyclic order of Z(VLIST)
  813.      remains unchanged, but the absolute order may change.
  814.  
  815.      This sorting order produces correct or nearly correct order
  816.      for a plfp command to make a plot involving hidden or partially
  817.      hidden surfaces in three dimensions.  It works best when the
  818.      polys form a set of disjoint closed, convex surfaces, and when
  819.      the surface normal changes only very little between neighboring
  820.      polys.  (If the latter condition holds, then even if sort3d
  821.      mis-orders two neighboring polys, their colors will be very
  822.      nearly the same, and the mistake won't be noticeable.)  A truly
  823.      correct 3D sorting routine is impossible, since there may be no
  824.      rendering order which produces correct surface hiding (some polys
  825.      may need to be split into pieces in order to do that).  There
  826.      are more nearly correct algorithms than this, but they are much
  827.      slower.
  828.  
  829.    SEE ALSO: get3_xy
  830.  */
  831. {
  832.   /* first compute z, the z-centroid of every poly
  833.    * get a list the same length as x, y, or z which is 1 for each
  834.    * vertex of poly 1, 2 for each vertex of poly2, etc.
  835.    * the goal is to make nlist with histogram(nlist)==npolys */
  836.   nlist= histogram(1+npolys(psum))(1:-1);
  837.   nlist(1)+= 1;  /* another problem with 1-origin indexing */
  838.   nlist= nlist(psum);
  839.   /* now sum the vertex values and divide by the number of vertices */
  840.   z= histogram(nlist, double(z))/npolys;
  841.  
  842.   /* sort the polygons from smallest z to largest z
  843.    * npolys(list) is the sorted list of lengths */
  844.   list= sort(z);
  845.  
  846.   /* next, find the list which sorts the polygon vertices
  847.    * first, find a list vlist such that sort(vlist) is above list */
  848.   vlist= list;
  849.   vlist(list)= indgen(numberof(list));
  850.   /* then reset the nlist values to that pre-sorted order, so that
  851.    * sort(nlist) will be the required vertex sorting list */
  852.   nlist= vlist(nlist);
  853.   /* the final hitch is to ensure that the vertices within each polygon
  854.    * remain in their initial order (sort scrambles equal values)
  855.    * since the vertices of a polygon can be cyclically permuted,
  856.    * it suffices to add a sawtooth function to a scaled nlist to
  857.    * produce a list in which each cluster of equal values will retain
  858.    * the same cyclic order after the sort
  859.    * (note that the more complicated msort routine would leave the
  860.    *  clusters without even a cyclic permutation, if that were
  861.    *  necessary) */
  862.   nmax= max(npolys);  /* this must never be so large that
  863.                * numberof(npolys)*nmax > 2e9  */
  864.   vlist= sort(nmax*nlist + indgen(numberof(nlist))%nmax);
  865.   /* primary sort key ^            secondary key  ^  */
  866. }
  867.  
  868. /* ------------------------------------------------------------------------ */
  869.  
  870. func spin3(nframes, axis, tlimit=, dtmin=, bracket_time=)
  871. /* DOCUMENT spin3
  872.          or spin3, nframes
  873.          or spin3, nframes, axis
  874.      Spin the current 3D display list about AXIS over NFRAMES.  Keywords
  875.      tlimit= the total time allowed for the movie in seconds (default 60),
  876.      dtmin= the minimum allowed interframe time in seconds (default 0.0),
  877.      bracket_time= (as for movie function in movie.i)
  878.  
  879.      The default AXIS is [-1,1,0] and the default NFRAMES is 30.
  880.  
  881.    SEE ALSO: rot3
  882.  */
  883. {
  884.   require, "movie.i";
  885.   if (is_void(nframes)) nframes= 30;
  886.   dtheta= 2*pi/(nframes-1);
  887.   if (is_void(axis)) axis= [-1.,1.,0.];
  888.   theta= acos(axis(3)/abs(axis(1),axis(2),axis(3)));
  889.   phi= atan(axis(2),axis(1)+(!axis(1)&&!axis(2)));
  890.   orig= save3();
  891.   movie, _spin3, tlimit, dtmin, bracket_time;
  892.   restore3, orig;
  893. }
  894.  
  895. func _spin3(i)
  896. {
  897.   if (i>=nframes) return 0;
  898.   rot3,,,-phi
  899.   rot3,,-theta,dtheta;
  900.   rot3,,theta,phi;
  901.   draw3;
  902.   return 1;
  903. }
  904.  
  905. /* ------------------------------------------------------------------------ */
  906.